Temporal and spatial correlations in a viscoelastic model of heterogeneous faults 
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We study the temporal and spatial correlations in a one-dimensional model of a heterogeneous 
fault zone, in the presence of viscoelastic effects. As a function of dynamical weakening and of 
dissipation, the system exhibits three different "phases" : one in which there are no time correlations 
between the events, a second, in which there are "Omori's law" type temporal correlations, and a 
third, runaway phase with quasiperiodic system size events. 
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I. INTRODUCTION 

The spatial and temporal distribution of earthquake 
activity has been the foremost aim and most successful 
aspect of earthquake modelling via discerctc, nonlinear 
networks of clastic elements interacting via nearest neigh- 
bor couplings, typically loaded at a constant rate far from 
the fault boundary. The Burridge-Knopoff 0] model has 
been the forerunner of a series of coarse-grained dynami- 
cal models , which have firmly established the under- 
standing of seismic activity within the paradigm of self- 
organized criticality These systems exhibit "sub- 
critical" or "supercritical" deviations J3| from strict self- 
similarity due to quenched inhomogeneities, finite driving 
velocities, or dissipation. 

A phenomenon which, to our knowledge, has not so 
far been addressed by dynamical models of the type 
cited above, is post-seismic relaxation M. It is com- 
monly believed that viscoelastic relaxation in the imme- 
diate post-seismic period results in a redistribution of the 
loads, with delay times of the order of minutes, hours, or 
days M. The modelling of these processes should help 
us understand such empirical findings as, for example, 
"Omori's Law," which says that the frequency of occur- 
rence of "aftershocks" decreases with time elapsed after 
the "main shock" as 



n(t) 



(const. + t)P 



(1) 



where p is usually found to be very close to unity |||To| . 

In this paper, we mimick viscoelastic relaxation by in- 
troducing a finite stress transfer velocity into a dynamical 
model recently studied by Dahmen et al. [fl3| , and we in- 
vestigate its effects on the spatio-temporal behaviour of 
this simple model. In our coarse grained representation, 
we do not claim to model the precise microscopic mech- 
anism for viscoelasticity, i.e., whether the relatively slow 
stress transfer actually comes from multiple brittle pro- 
cesses in a heterogenous medium | pT[ or from coupling to 
a viscous layer below the litosphere Jl2j . We will simply 
take the stress transfer velocity (V) to be some effec- 
tive group velocity which governs post-event relaxation 



in the system and which is smaller than the velocity of 
sound @. 

The model we have used as our point of departure is 
an infinitely long range (Mean Field) version of the Ben- 
Zion and Rice |I| model which has been investigated both 
analytically and numerically |l5| ], to reveal the presence 
of two different regimes as far as spatial and temporal 
distributions are concerned. It has been found, for a 
narrow distribution of heterogeneities, in the limit of in- 
finitely slow drive, that the phase space can be described 
in terms of just two parameters, the dynamical weaken- 
ing e and conservation c, both taking values between 
and 1. For c < 1/2 and small e, the behaviour is critical; 
this is the so called Gutenberg-Richter (GR) regime, with 
a power law distribution of event sizes. For c > 1/2 and 
e close to unity, one finds a metastable state of two-phase 
coexistance, with GR behaviour interrupted by stretches 
of quasi-periodic, characteristic (system-size) events, i.e. 
"runaway" behaviour. 

Clearly, for the purely Abelian models that have so far 
been considered |]7|,|l3|, just the retardation effect com- 
ing from the introduction of a finite velocity of stress 
transfer cannot make any difference in the overall statis- 
tics, since it does not matter in which sequence the sites 
are updated jn|. Therefore, the introduction of spatial- 
ity beyond that embodied in the mean field (infinitely 
long range) approximation had to be considered simul- 
taneously with the retardation effect. Namely, the in- 
teractions strengths ("spring constants") were made to 
depend inversely on the distance, in keeping with a one- 
dimensional picture of the fault zone. 

In this study we therefore consider simultaneously 
zjthe effect of time delay in the transfer of stress, ii) the 
decay of the coupling strength with inverse distance, for 
the model of Dahmen et al. [H| in one dimension. We 
have moreover considered a slightly different distribution 
of heterogeneities. 

We find that this new model leads to three distinct 
phases. The phase diagram is shown schematically in 
Fig. 1. The phase space has been probed over a grid of 
Ac = 0.1, for e = 0, 0.5 and 1. For strongly dissipative (c 



1 



near 0) systems with relatively weak "dynamical soften- 
ing" effects (e close to 0) we find a GR-like phase, with 
very small events which show a very steep incipient power 
law behaviour over a rather narrow range of sizes and 
then cut off abruptly, and which display essentially no 
temporal correlations. The power spectrum of the event 
sequence is white-noise. For intermediate values of these 
parameters, the event distribution is similar to the first, 
however the power spectrum reveals non-trivial temporal 
correlations, a feature not observed by Dahmen et al. |f| 
in the GR phase. In the region of c close to unity (strong 
conservation)we find quasiperiodic runaway behaviour. 

1-S 




FIG. 1. Schematic phase diagram exhibiting three differ- 
ent "phases" where, I.) small events delta correlated in time, 
II.) small events with Cauchy type time correlations, and III.) 
quasiperiodic, system-size runaway events, dominate. 

The phenomenology of each of these phases is rather 
rich. In regions of strong heterogeneity, one may ob- 
serve patches of blocks to slip in unison, pinned at ei- 
ther end by relatively large threshold stresses, exhibiting 
quasiperiodic behaviour within a sea of power law events. 
We have not observed the switching behaviour between 
coexisting GR and runaway phases, as reported for the 
mean field model Jl3| , but this may be because this would 
require prohibitively long simulations in our scheme. The 
distribution of the accumulated stresses along the fault 
zone is both qualitatively and quantitatively similar in 
all the three regions, in contrast to the previous findings. 

It should be stressed here that our approach is a de- 
parture from the usual quest for scale invariant spatial 
or temporal distributions. With the introduction of a fi- 
nite stress transfer velocity (which we take to be unity), 
and a finite driving velocity, we in fact have three well 
seperated time scales in the problem: That of the driv- 
ing velocity (the largest time scale) , the viscoelastic time 
scale, and the triggering time scale (where slip occurs 
instantaneously) . 

The paper is organized as follows. In section 2, the 



precise definiton of the model is given. In section 3, we 
report our results for the statistics of the magnitudes in- 
tegrated over time scales corresponding to typical event 
durations, which should be compared with those in the 
infinite stress transfer velocity/zero driving velocity limit. 
We then go on to compute temporal and spatial correla- 
tion functions for coarse grained events. In section 4 we 
provide a discussion of our findings. 



II. THE MODEL 

We consider a one-dimensional array of finite segments, 
or blocks. The local stress r, on the ith block is given, 
at time t by 

R 

n{t)= k r [u i+ r(t-r/V)-Ui(t)]+K[vt-Ui(t)] 



r=-R 



(2) 



where the range of the interaction, R, is of the order of 
the system size, Ui{t) is the offset of the zth block in the 
direction of the constant driving velocity v, at time t; 
K is the effective shear modulus, and k r — k/\r\ is the 
elastic coupling between blocks seperated by a distance 
r. As long as all the Tj < r St i 1 where {t s ^} are randomly 
distributed failure stresses, the system is immobile. 

The viscoelastic stress relaxation is mimicked by the 
delay, r/V, in the transfer of stress. We shall henceforth 
set V, the velocity for the stress transfer along the blocks, 
to unity. Note that V is not typically the sound veloc- 
ity, but some effective group velocity smaller than that 
of sound, governing the processes of viscoelastic stress 
relaxation in this coarse grained model. |8HI^] 

The dynamics is defined as follows. If the threshold 
value is exceeded at some i, at time t, then 

i)The stress at the ith block is reduced by 



5n 



(3) 



where the {r a ,i} are random arrest stresses. 

ii) The value of the failure stress at the i'th block is re- 
set, until all motion once more ceases, to a "dynamical" 
threshold value 



T,l. 



^(Ts,2 7~a,i ) ) 



(4) 



where e parameterizes the dynamical weakening effect. 

in) The stress drop is redistributed, according to 
Eq.(|J), so that r^ +r is incremented, at the t+ |r|'th time 
step by 



K 



(5) 



We may once more define c = ^ r c r , with < c < 1, to 
be the parameter which measures the degree of conserva- 
tiveness of the system, although it should be noted that 



2 



this definition now involves an implicit integral over time 
as well as space. 

The boundary conditions are fixed, so that u% = Ul = 
0. At each time step, the stress at all the blocks i is 
recalculated according to Eq.(||). This means that the 
constant drive term Kvt is incremented also. 



III. SIMULATIONS 

Since the finite stress transfer velocity introduces a def- 
inite time scale into the system, which also sets the char- 
acteristic time scale of the event duration, one now faces 
the problem of having to go to extremely long runs with 
a driving velocity which is at least six to seven orders 
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FIG. 2. Typical time series for the magnitudes M 
v.s. time, for 16384x128 time steps, a) in the "Guten- 
berg- Richter" (small event) and b) the runaway regimes. Note 
the difference in the vertical scales. Shown are plots for a) 
c = 0.16, e = 0.5; b) c = 0.9, e = 0.5. The driving velocity 
is v — 10~ 5 and the shear modulus K = 1 for this and the 
following figures. 

of magnitude smaller than the latter. The zero driving 
velocity trick of simply scanning the system for that site 
which is closest to slipping and loading all sites by the 
missing amount, is no longer appropriate here - one has 
to first check that there are no stress "parcels" still on 



the way. The actual simulation times get prohibitively 
large as a result, and we had to be content with large 
but finite ranges of interaction, up to 1/6 the fault size, 
and with a one-dimensional fault. 

Table I. 

Values of the Gutenberg-Richter exponent b in the dif- 
ferent regimes. 



Region 


c 


e 


b 


I 


0.15 


0.5 


5.2 


II 


0.65 


0.5 


5.2 


III 


0.9 


0.5 


1.6 



We have simulated the system described in the pre- 
ceding section on a grid of 300 blocks, with the range 
of interactions going up to R = 50. The distribution of 
stress drops 5ti = t s j — r a ^ was chosen to have the form 
p(x) oc x~^j with /j, = 1.2 However, upon finding that ar- 
bitrarily large stress drops pinned the edges of finite seg- 
ments in the fault and distorted the distribution of event 
sizes, we decided to limit the range of the St to a width 
comparable to those considered in Ref. ( jlH), namely 
0.2. It is generally found that long active fault zones 
organize themselves into states with relatively small het- 
erogeneity, and our results should be considered in this 
spirit. 

The driving velocity v is taken to be 10~ 5 V in the sim- 
ulations reported below. It should be noted that larger 
driving velocities result in individual cells exceeding their 
threshold values and collapsing independently from their 
neighbors, and as a result the system never achieving a 
self organized state. With realistic driving velocities of 
~ 10~ 9 m/sec, v = 10~ 5 V corresponds to a stress trans- 
fer velocity of ~ 10~ 4 m/sec. For "block sizes" of ~ 10 2 
meters setting our lattice spacing, a stress transfer veloc- 
ity of V ~ 10 m/s corresponds to time steps of duration 
10 6 s. 

In Figure 2, we display the typical time series resulting 
from plotting the integrated magnitudes M(t) for differ- 
ent values of the system parameters. These are obtained 
by summing over the number of blocks where the thresh- 
old has been exceeded within an interval St — 128. 
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FIG. 3. The frequency v.s. magnitude plots on a double 
logarithmic scale, for the same set of parameters as a) and b) 
in Figure 2. The lines are intended as a guide to the eye with 
slopes —5.2 and —1.6, respectively. 

This value was chosen as approximately that time in- 
terval needed for a signal originating in the middle of the 
fault to be able to reach the edges of the fault. This 
time-coarse- grained way of identifying events is in keep- 
ing with the way earthquake data is taken, with the time 
integral taken over the actual displacements of the seis- 
mographs. Thus, we define, 



t+st 

M(t) = £5} Au *(*')]° 



(6) 



where the zeroth power of the slip Au,(i') is taken so 
that M(t) simply counts the number of slipped blocks 
within the time interval St. The time series M{t) re- 
veal no immediately observable differences between the 
regions I and II shown in Fig.l; therefore we have se- 
lected only one set of parameter values to illustrate both 
these regions (see Fig. 2a). On the other hand the much 
larger magnitudes observed in the runaway region, and 
their marked quasiperiodicity are apparent in Fig. 2b. 

The regions I and II are also similar in the way the fre- 
quency f(M) of events scales with the magnitude M, for 



a given binning size St. In Figure 3, we show the plots of 
the frequency f(M) v.s. magnitude M in the small event 
(the "Gutenberg-Richter" phase found in Ref. and 
in the runaway regime, for the same parameter values 
as in Fig. 2. The "power law" fits, f(M) ~ M~ b , to 
the small event regime (regions I and II) are poor, and 
can only be thought of as suggestive; they extend over 
too small a range to really signify self-similarity. Notice 
that in the "runaway" regime (region III), the magni- 
tudes cover a wider range; however we do not observe 
as marked a "super-criticality," i.e., a frequency of large 
events in excess of a power law size distribution, as has 
been reported elsewhere p,0JlT 
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FIG. 4. Power spectra of the time series of magnitudes, 
computed for parameter values in the regions I, II, III of the 
phase diagram shown in Fig.(l). a) c = 0.16, e = 0.5, b) 
c = 0.72, e = 0.5, c) c = 0.89, e = 0.5. The data for panels 
a) and c) are taken over series of 8192x128 time steps. The 
data in panel (b) have been averaged over 7 runs of 16384x 128 
time steps. 

Within regions I and II, the Gutenberg-Richter expo- 
nent b is sensitive to the driving velocity v and also de- 
pends, less strongly, on the parameters c, e. Here the dis- 
tribution is very steep, with b ranging between 4 and 5. 
The values obtained are given in Table I. The "runaway" 
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regime exhibits much more realistic b values, around 1.5 
to 1.6. One should note, moreover, that since the system 
is no longer scale invariant, the statistics of the magni- 
tudes M are also sensitive to the binning size St, so that 
the b values here are only useful for purposes of discrim- 
inating between different regions of the phase space. 

The interesting difference between the three regimes 
delienated in Fig. 1 become apparent in the power spec- 
tra, 



S(f) = J dte a ^C(t) , 



(7) 



where 



C(t) = 



dt'M(t')M(t' + 1) 



(8) 



is the time-correlation function for the coarse grained 
magnitudes M(t), where t stands for the number of time 
intervals St, for a total time of measurement extending 
over a period T. 
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FIG. 5. The same power spectra as in Fig. 4, on a log- 
arithmic scale. Panels a), b), and c) belong to the regions I, 
II, III of the phase diagram in Fig. 1. 



In Figures 4 and 5 we display linear and semi- 
logarithmic plots of the power spectrum in the three dif- 
ferent phases. It can be seen in Fig. (4a) and in a more 
pronounced way in the logarithmic plot of the same power 
spectrum in Fig. (5a) that in the "small event" regime I 
(panel a), the power spectrum is essentially flat, white- 
noise like, indicating an absence of correlations between 
the earthquakes, i.e., C(t) oc S(t). For intermediate val- 
ues of c and e, i.e., in the small-event region II, how- 
ever, we find that superposed upon the white-noise like 
background, the upper envelope of the power spectrum 
displays a distinctive curve, (see Fig. 4b) indicating the 
presence of non-trivial temporal correlations. In the run- 
away region (region III), we find a markedly different, 
quasiperiodic behaviour, as can be see from Figs. (4c) 
and (5c). The very pronounced peak in the power spec- 
trum near the origin is large enough to suppress all the 
others; we can see the other frequencies that are present 
only in the logarithmic plot. 

In Figure 6, we show the result of taking an inverse 
transform of the envelope (roughly the highest points) of 
the power spectrum shown in Fig. 4b. This crude esti- 
mate of the time correlation function is corraborated by 
a more careful evaluation, to which we now turn. 




FIG. 6. A crude estimate for the time correlation function, 
from the inverse Fourier transform of the points in the enve- 
lope of the power spectrum, shown in Fig. 4b. The fit is to 
1 + 85(1 + St)- 1 . 

We have computed C(t) directly from a time series of 
81920 time intervals (of 128 steps each). We have nor- 
malized the correlation function by (1/T) M 2 (t')dt' , 
so that C(0) = 1. We find that the normalised C(t) can 
be fit rather well by a function of the form 



C(t) = A + 



B 



D + t 



(9) 
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FIG. 7. The time correlation function C(t) in region II 
of the phase diagram (c = 0.65, e = 0.5), computed from 
Eq.(8), over a time series of 81920 x 128 steps. The fit is to 
0.7 + 0.2(0.66+ t)' 1 . 
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FIG. 8. The probability density of the fraction s of the 
slipping stress (see text) computed along the fault zone, for 
c = 0.89, e = 0.5 (region III). The histogram is averaged over 
20 snapshots, seperated by 5 x 10 4 timesteps. The figures for 
regions I and II are indistinguishable from this one. 

To investigate the spatial correlations in the system, 
we first considered the distribution of the fraction of the 
slipping stress on each block, 



with A — 0.7, B — 1/5, D — 2/3. Our results are shown 
in Fig. 7. 

Note that the time correlation function C(t) measures 
the average frequency with which a time lapse t seperates 
two events, weighted by the magnitude of the events. In 
another way of saying the same thing, it is the weighted 
average of the number of times that one registers pairs 
of events seperated by a time lapse equal to t. If there 
is no event taking place at time t + if after an event at 
time i.e., if M(t') ^ 0, but M(t + t') = 0, there will 
be no contribution to the integral for C(t) in (|§|). 

We would like to recall, at this point, Omori's Law 
(|I|) for the frequency of aftershocks |3 10|. Actually, geo- 
physicists are generally hesitant to label a given shock as 
either an "aftershock" or "main schock," and admit that 
these are conventional distinctions, which are difficult to 
make precise in a quantitative way. Viewed in this way, 
Omori's law is just a statement of the relative frequency 
of pairs of events seperated by a time t, and is a slightly 
cruder version of the time correlation function. We see 
that the form we find for the time-decay of the correla- 
tion function matches that of Omori's Law, with a power 
p = 1, as found for real earthquake statistics. 

From Fig. 7, and Eq.(9), we see that C(t) ~ A drops 
by a factor of 1/2 within one time interval (consisting 
of 128 time steps). Since we have already estimated the 
time steps here to correspond to about 10 6 seconds, this 
means correlations times of the order of 10 8 s, namely 
~ 3 yrs, which is quite realistic for the time interval in 
which aftershocks die away after a big event. 



(10) 



along the fault zone. This is a quantity in which one 
might have expected a greater self-organization building 
up as one goes from region I to region III in the phase 
diagram, yet we did not find this to be the case. In con- 
trast to Ref. 13 1, in the present model the distribution 
p(s) remains essentially invariant in all the three regions, 
and has the shape shown in Fig. 8. Similarly, wc found 
that the equal time spatial correlations {siSi+ r } between 
the fraction of the slipping stresses accumulated at each 
site, showed essentially delta function behaviour in all 
three regions, with the correlations never extending be- 
yond next nearest neighbors. On the other hand, defining 
the coarse grained toppling variables 



t+st 
t'=t 



(11) 



and setting m c (i) — at sites beyond the boundaries of 
the fault we found that the correlation function 



C m (r) 



(m c (t, i)m c (t, i + r)) - (m c (t, i)) 2 
(m c (t, i)) 2 



(12) 



where the averages are performed both over i and t, in- 
deed displayed markedly different behaviour in region III, 
in comparison to I and II. Our results are shown in Fig. 
9. The time averages were performed over 60 snapshots, 
taken at intervals of 5 x 10 4 timesteps. The correlations 
are negligible (of the order of 10~ 4 ) in the first two re- 
gions, whereas, in region III, one sees a gradual decay. A 
straight line fit to the semilogarithmic plot (Fig. 10) sug- 
gests an exponential decay and gives a correlation length 
of 30 lattice units, corresponding to ~ 3 x 10 3 meters. 
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FIG. 9. The spatial correlations between toppling events 
coarse grained in time, in the regions I and III of the phase 
diagram, a) c = 0.16, e = 0.5, b) c = 0.89, e = 0.5. The 
plot for c = 0.54,e = 0.5 (in region II) looks identical to panel 
a). Averages have been taken over 60 snapshots seperated by 
time intervals of 50,000 steps. 



IV. DISCUSSION 

The inclusion of viscoleastic effects into the study of 
crack propagation and pinned driven systems has re- 
cently made important progress [ tLaj^ , and promises 
to be fruitful also in the modeling of earthquakes. The 
most important outcome of introducing viscoelastic ef- 
fects is to be found in the more subtle temporal correla- 
tions between events; for highly dissipative systems the 
correlations are delta function like, whereas for interme- 
diate values of the dissipation, one observes correlations 
that decay as ~ (const. + t)^ 1 between events, which is 
of the form of Omori's Law ([!]). To our knowledge, this 
is the first demonstration of how Omori's Law may arise 
in such a system. 

Our preliminary findings indicate that, due to vis- 
coelastic effects, the runaway phase (region III) of 
quasiperiodic events in the present model of a heteroge- 
neous fault zone is pushed to a relatively smaller region of 



the phase diagram (see Fig. 1) than found previously (51| 
and the frequency distribution in this phase displays scale 
invariance over a sizable region of event sizes. We have 
verified that this region is distinguished by relatively long 
range spatial correlations between slipping events, in con- 
trast to the "small event" regimes. 



FIG. 10. Semilogarithmic plot of the spatial correlations 
in region III. The slope of the straight line fit gives a correla- 
tion length of 30 lattice units, or equivalently, 3 x 10 meters. 

It has been remarked before p,[I8| that various system- 
dependent features, notably dynamical weakening and 
dissipation, introduce time and length scales into the 
problem and take the system away from criticality. We 
would like to remark that one may reverse the empha- 
sis here to say that rather than the scaling region of the 
Gutenberg-Richter regime, one should examine the sub- 
or super-critical behaviour to characterise a specific fault 
zone. In particular, we have shown that determining the 
nature of the space and time-correlations in the system 
gives important clues as to the relative degree of dissipa- 
tion or dynamical weakening. 
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